library(Seurat)
library(tidyseurat)
library(harmony)
library(SingleR)
library(ggpubr)
library(ggrepel)
library(tidyverse)
source('00_util_scripts/mod_seurat.R')
source('00_util_scripts/mod_bplot.R')

sobj <- read_rds('mission/ye_AH-cOME/ah9.rds')

# remove T & mDC
sobj <- sobj |>
  filter(!(RNA_snn_res.0.8 %in% c(12, 14)))

sobj <- sobj |>
  quick_process_seurat(batch = c('orig.ident','sex'), pcs = 50,
                       skip_norm = TRUE)

sobj |>
  DotPlot(c('S.Score', 'G2M.Score'))

# annotate with markers in Wang-Ye paper
ye.bcell.type <- list(
  naive.B = c('TCL1A','SELL','FCER2'),
  FCRL4.MBC = c('SOX5','FCRL4','FCRL5','PTPN1','PLAC8','ITGAX'),
  classic.MBC = c('TNFRSF13B','BANK1','CD27'),
  LZ.GC = c('CD83','BCL2A1'),
  DZ.GC = c('CD38','SUGCT','AICDA'),
  cycling = c('HMGB2','MKI67','TUBA1B'),
  plasmablast = c('XBP1','MZB1','JCHAIN'),
  CD21low = c('COX5A','MYL6','UQCRH','COTL1')
)

sobj |>
  DotPlot(ye.bcell.type, cluster.idents = T) +
  rotate_x_text()

sobj <-
"
RNA_snn_res.0.8,manual_wang
1,Naive.B
2,DZ.GC
3,cycling
4,cycling
5,FCRL4.MBC
6,plasmablast
7,CD21low
8,LZ.GC
9,cycling
10,DZ.GC
11,classical.MBC
12,plasmablast
" |> read_delim() |>
  mutate(RNA_snn_res.0.8 = as_factor(RNA_snn_res.0.8)) |>
  left_join(sobj, y = _)

sobj |>
  DimPlot(group.by = 'manual_wang', cols = 'Set3', label = T) +
  ggtitle('Adenoid B cells subsets')

sobj |>
  DotPlot(list_c(ye.bcell.type), group.by = 'manual_wang', cluster.idents = T) +
  RotatedAxis()

sobj |>
  mutate(cOME = ifelse(cOME == 1, 'cOME', 'Ctrl')) |>
  write_rds('mission/ye_AH-cOME/ah9.rds')

# annotate with ref data ------
sobj |> DimPlot(group.by = 'hpca.fine')

sobj |> FeaturePlot(c('FCRL4','AICDA','CD38'))

latent.come <- sobj |>
  annotate_latents('B.cell', logfc.thres = 1.5)

# compare cell fraction between group -----------
## cOME vs ctrl --------
frac.come <- sobj@meta.data |>
  discov_frac_change(cOME, manual_wang, cOME, Ctrl)

frac.come |>
  ggplot(aes(subtype, log2fc_frac, fill = type)) +
  geom_col() +
  scale_fill_manual(values = c('blue','red','grey')) +
  theme_pubr() +
  coord_flip() +
  labs(title = 'B cell subset fraction changes in AH + cOME vs AH patients')

sobj@meta.data |>
  dplyr::count(orig.ident, manual_wang, cOME) |>
  group_by(orig.ident) |>
  reframe(cOME ,manual_wang, fraction = n / sum(n)) |>
  mutate(subtype = manual_wang, .keep = 'unused') |>
  left_join(frac.come[,c(1,5)]) |>
  mutate(subtype = fct_reorder(subtype, log2fc_frac, .desc = T),
         cOME = ifelse(cOME == 'Ctrl', 'AH','AH+cOME') |> fct_relevel('AH')) |>
  ggplot(aes(cOME, fraction, fill = cOME)) +
  stat_summary(geom = 'col') +
  geom_jitter(width = .01, height = 0) +
  facet_wrap(~subtype, scales = 'free_y') +
  theme_pubr() +
  scale_fill_manual(values = c('skyblue','red')) +
  labs(title = 'B cell subset fraction changes by sample',
       fill = 'group', x = 'group')

## RR vs GG+GR --------
frac.rr <- sobj@meta.data |>
  mutate(RRvG = ifelse(`IGHG1-G396R` == 'RR', 'RR', 'other')) |>
  discov_frac_change(RRvG, manual_wang, RR, other)

frac.rr |>
  ggplot(aes(subtype, log2fc_frac, fill = type)) +
  geom_col() +
  scale_fill_manual(values = c('blue','red','grey')) +
  theme_pubr() +
  coord_flip() +
  labs(title = 'B cell subset fraction changes in RR vs GG+GR patients')

## IT vs II -----
frac.it <- sobj@meta.data |>
  discov_frac_change(`FCGR2B-I232T`, manual_wang, IT, II)

frac.it |>
  ggplot(aes(subtype, log2fc_frac, fill = type)) +
  geom_col() +
  scale_fill_manual(values = c('blue','red','grey')) +
  theme_pubr() +
  coord_flip() +
  labs(title = 'B cell subset fraction changes in IT vs II patients')

# DEG ---------
## between AH+COME vs AH --------
Idents(sobj) <- 'manual_wang'

come.type.list <- unique(Idents(sobj))

come.type.list <- come.type.list |>
  set_names(come.type.list)

deg.subsets <- come.type.list |>
  map(\(x)FindMarkers(sobj, ident.1 = 'cOME', group.by = 'cOME',
                      subset.ident = x) |>
        as_tibble(rownames = 'gene'), .progress = T) |>
  list_rbind(names_to = 'subtype')

sdeg.subsets <- deg.subsets |>
  select(-p_val) |>
  filter(p_val_adj < .05) |>
  write_csv('mission/ye_AH-cOME/deg.subsets.csv')

sdeg.subsets <- sdeg.subsets |>
  mutate(p_val_adj = ifelse(p_val_adj == 0, 1e-300, p_val_adj))

volc.list <- come.type.list |>
  map(\(x)filter(sdeg.subsets, subtype == x, str_starts(gene, 'MTND|ENSG|LINC', negate = T)) |>
        plot_bill_volc('AH+cOME') +
        ggtitle(x))
  
volc.list[[8]]

## between G2R --------
### RR vs GG+GR
deg.subsets.g2r <- come.type.list |>
  set_names(come.type.list) |>
  map(\(x)FindMarkers(sobj, ident.1 = 'RR', group.by = 'IGHG1-G396R',
                      subset.ident = x) |>
        as_tibble(rownames = 'gene'), .progress = T) |>
  list_rbind(names_to = 'subtype')

sdeg.subsets.g2r <- deg.subsets.g2r |>
  select(-p_val) |>
  filter(p_val_adj < .05) |>
  write_csv('mission/ye_AH-cOME/deg.subsets.g2r.csv')

### GG vs GR+RR
deg.subsets.gg <- come.type.list |>
  set_names(come.type.list) |>
  map(\(x)FindMarkers(sobj, ident.1 = 'GG', group.by = 'IGHG1-G396R',
                      subset.ident = x) |>
        as_tibble(rownames = 'gene'), .progress = T) |>
  list_rbind(names_to = 'subtype')

sdeg.subsets.gg <- deg.subsets.gg |>
  select(-p_val) |>
  filter(p_val_adj < .05) |>
  write_csv('mission/ye_AH-cOME/deg.subsets.gg.csv')

## RR vs GG
deg.subsets.rrgg <- come.type.list |>
  set_names(come.type.list) |>
  map(\(x)FindMarkers(sobj, ident.1 = 'RR', ident.2 = 'GG', group.by = 'IGHG1-G396R',
                      subset.ident = x) |>
        as_tibble(rownames = 'gene'), .progress = T) |>
  list_rbind(names_to = 'subtype')

sdeg.subsets.rrgg <- deg.subsets.rrgg |>
  select(-p_val) |>
  filter(p_val_adj < .05) |>
  write_csv('mission/ye_AH-cOME/deg.subsets.rrgg.csv')

## between I2T --------
deg.subsets.i2t <- come.type.list |>
  set_names(come.type.list) |>
  map(\(x)FindMarkers(sobj, ident.1 = 'IT', group.by = 'FCGR2B-I232T',
                      subset.ident = x) |>
        as_tibble(rownames = 'gene'), .progress = T) |>
  list_rbind(names_to = 'subtype')

sdeg.subsets.i2t <- deg.subsets.i2t |>
  select(-p_val) |>
  filter(p_val_adj < .05) |>
  write_csv('mission/ye_AH-cOME/deg.subsets.i2t.csv')


# try assign IG isotype B cell
ig.isotype <- sdeg.subsets |>
  filter(str_detect(gene, '^IGH[D|M|A|E|G]')) |>
  pull(gene) |> unique()

cb.isotype <- sobj |>
  get_abundance_sc_long(ig.isotype, exclude_zeros = T)

pure.isotype <- cb.isotype |>
  dplyr::count(.cell) |>
  filter(n == 1) |>
  left_join(cb.isotype) |>
  mutate(isotype = .feature, .cell = .cell, .keep = 'none')
  
pure.isotype <- sobj@meta.data |>
  as_tibble(rownames = '.cell') |>
  right_join(pure.isotype)

pure.isotype |>
  mutate(cOME = ifelse(cOME == 'cOME', 'AH+cOME','AH')) |>
  ggplot(aes(cOME, fill = isotype)) +
  geom_bar(position = 'fill') +
  coord_flip() +
  scale_fill_brewer(palette = 'Paired') +
  theme_pubr() +
  labs(title = 'B cell isotype fraction of AH+cOME and AH patients',
       y = 'Fraction',
       x = 'Group')

pure.isotype |>
  ggplot(aes(`IGHG1-G396R`, fill = isotype)) +
  geom_bar(position = 'fill') +
  coord_flip() +
  scale_fill_brewer(palette = 'Paired') +
  theme_pubr() +
  labs(title = 'B cell isotype fraction of AH patients with IGHG1-G396R genotype',
       y = 'Fraction',
       x = 'IGHG1-G396R genotype')

pure.isotype |>
  ggplot(aes(`FCGR2B-I232T`, fill = isotype)) +
  geom_bar(position = 'fill') +
  coord_flip() +
  scale_fill_brewer(palette = 'Paired') +
  theme_pubr() +
  labs(title = 'B cell isotype fraction of AH patients with FCGR2B-I232T genotype',
       y = 'Fraction',
       x = 'FCGR2B-I232T genotype')

# key cytokine violin --------
sobj |> VlnPlot(c('TNF','IL1B','IL4','IL6','CXCL8','IL17A'),
                pt.size = 0)

sobj |>
  filter(str_detect(manual_wang, 'MBC')) |>
  VlnPlot('TNF', group.by = 'IGHG1-G396R') +
  labs(title = 'TNF expression in memory B cells from AH patients',
       fill = 'IGHG1-G396R genotype', x = 'IGHG1-G396R genotype')

sobj |>
  filter(str_detect(manual_wang, 'MBC')) |>
  mutate(cOME = ifelse(cOME == 'cOME','AH+cOME','AH')) |>
  VlnPlot('TNF', group.by = 'cOME') +
  labs(title = 'TNF expression in memory B cells from AH patients',
       x = 'Group')

# allele exlcuded G2R B cell in GR individual -------
bc.g2r.call <- read_csv('mission/ye_AH-cOME/ah-come.g2r.bc.csv')

bc.g2r.call <- bc.g2r.call |>
  mutate(.cell = str_c(pid, '_', .cell))

bc.g2r.call |>
  summarise(n = n(), .by = .cell) |>
  filter(n > 1) |>
  left_join(bc.g2r.call) |>
  slice_max(baseq, by = .cell, with_ties = F) |>
  mutate(g2r.AA = case_match(g2r.base, 'C' ~ 'only G',
                               'T' ~ 'only R', .default = 'Both G&R')) |>
  ggplot(aes(g2r.AA, fill = g2r.AA)) + geom_bar()

al.ex.gr <- bc.g2r.call |>
  filter(str_detect(g2r.base, 'X', negate = T), baseq > 10) |>
  left_join(ah.meta) |>
  mutate(.cell = str_c(id, '_', .cell), g2r.base, mapq, baseq,
         .keep = 'none') |>
  slice_max(baseq, by = .cell, with_ties = F) |>
  left_join(x = sobj, y = _) |>
  filter(!is.na(g2r.base))

al.ex.gr |>
  mutate(g2r.aa = ifelse(g2r.base == 'C', 'G', 'R')) |>
  ggplot(aes(g2r.aa, fill = manual_wang)) +
  geom_bar(position = 'fill') +
  scale_fill_brewer(palette = 'Paired') +
  coord_flip() +
  labs(y = 'Fraction', x = 'AA residue on IGHG1-G396R', fill = 'B cell subtype')

sdeg.al.ex.g2r <- al.ex.gr |>
  FindMarkers(group.by = 'g2r.base', ident.1 = 'T') |>
  as_tibble(rownames = 'gene')

sdeg.al.ex.g2r |>
  filter(p_val_adj < .05)

sdeg.al.ex.g2r |>
  mutate(p_val_adj = p_val) |>
  plot_bill_volc(exp_group = 'R') +
  ylab('-log10(p)')
